Multiscale modelling of liquids with molecular specificity 
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The separation between molecular and mesoscopic length and time scales poses a severe limit 
to molecular simulations of mesoscale phenomena. We describe a hybrid multiscale computational 
technique which address this problem by keeping the full molecular nature of the system where it is 
of interest and coarse-graining it elsewhere. This is made possible by coupling molecular dynamics 
with a mesoscopic description of realistic liquids based on Landau's fluctuating hydrodynamics. We 
show that our scheme correctly couples hydrodynamics and that fluctuations, at both the molecular 
and continuum levels, are thermodynamically consistent. Hybrid simulations of sound waves in bulk 
water and reflected by a lipid monolayer are presented as illustrations of the scheme. 
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Complex multiscale phenomena are ubiquitous in na- 
ture in solid (fracture propagation Q), gas (Knudsen lay- 
ers 0) or in liquid phases (fluid slippage past surfaces 
, crystal growth from fluid phase, wetting, membrane- 
fluid dynamics, vibrational properties of proteins in wa- 
ter 0,0 and so on). These phenomena are driven by 
atomistic forces but manifest themselves at larger, meso- 
scopic and macroscopic scales which cannot be resolved 
by purely large scale molecular simulations (with some 
notable exceptions 0]). On the other hand, coarse- 
grained mesoscopic models have limited use due to the 
approximations necessary to treat the molecular scales 
intrinsic to these methods. A viable solution to this 
dilemma is represented by multiscale modelling via cou- 
pled models, a protocol which is also well suited to new 
distributed computing paradigms such as Grids 0, HI- 
The idea behind this approach is simple: concurrent, cou- 
pled use of different physical descriptions. 

The coupled paradigm is the underlying concept in 
quantum-classical mechanics hybrid schemes |l| used to 
describe fracture propagation in brittle materials and also 
in hybrid models of gas flow Q • During the last decade, 
hybrid modelling of liquids has received important con- 
tributions from several research groups (see the recent 
review |9|). However, it has thus far lacked the maturity 
to become a standard research tool for liquid and soft 
condensed matter systems. Hybrid simulations of liq- 
uids have been restricted to coarse-grained descriptions 
based on Lcnnard- Jones particles, reducing the major 
advantage of this technique of maintaining full molec- 
ular specificity where needed. Recently, new methods for 
energy controlled insertion of water molecules [Io| have 
finally opened the way to real solvents such as water. So 
far, no hybrid method has employed an accurate descrip- 
tion of the mesoscale (from nanometres to micrometres) 
as the important contribution of fluctuations has been 
neglected in the embedding coarse-grained liquid. The 
hybrid method must also ensure thermodynamic consis- 
tency, by allowing the open molecular system to relax to 
an equilibrium state consistent with the grand canonical 



ensemble 1 1 11 - Finally, all previous non-equilibrium hy- 
brid simulations have been restricted to shear flow |°Ill2fl. 

In this Letter, we present a coupled multiscale model 
called "hybrid MD" for simulation of mesoscopic quanti- 
ties of liquids (water) embedding a nanoscopic molecular 
domain (Fig. la). Hybrid MD overcomes the limitations 
of previous hybrid descri ptio ns of liquids by coupling fluc- 
tuating hydrodynamics |l3| and classical molecular dy- 
namics via a protocol which guarantees mass and momen- 
tum conservation. The present method is designed to ad- 
dress phenomena driven by interplay between the solute- 
solvent molecular interaction and the hydrodynamic flow 
of the solvent. 

Fluctuating hydrodynamics model. Our mesoscopic de- 
scription of fluid flow is based on the equations of fluc- 
tuating hydrodynamics (FH) ^jjj]. These equations are 
stochastic partial differential equations which reduce to 
the Navier-Stokes equations in the limit of large volumes. 
The equations are based on the conservation equations 
dt<t> — — VJ^, where c/> = 4'{ r ^) is the density of any 
conserved variable at location r. We consider an isother- 
mal fluid, so that the relevant variables are the mass and 
momentum densities 4> = {p, g} (here g = pv and v 
is the fluid velocity). The mass and momentum fluxes 
are given by J p = pv and J 9 = pvv + II + II, where 
II and II are the mean and fluctuating contributions 
to the pressure tensor, respectively. The mean pressure 
tensor is usually decomposed as II = (p + ir)l + II s , 
where p is the thermodynamic pressure (given by the 
equation of state) and the stress tensor is the sum of a 
traceless symmetric tensor II S and an isotropic stress 
7r. We consider a Newtonian fluid for which II = 
—r\ [poVp 4- dpv a — 2D~ 1 d 1 v~ l 8 a p) , ir = — £<5 7 v 7 , where 
repeated indices are summed, D the spatial dimension 
and 77, C are the shear and bulk viscosities respectively. 
The components of the fluctuating pressure tensor n Q( g 
are random Gaussian numbers (see supplementary infor- 
mation) . 

Our continuum mesoscopic model is based on a finite 
volume discretization of the FH equations 0] , although 



here in an Eulerian frame of reference and on a regular 
lattice. Partitioning the space into several space-filling 
volumes V k with k = 1, N centered at positions r k , we 
integrate the conservation equations over each volume V k 

and apply Gauss' theorem 4? J v <f>(rk,t)dr — Au^ti ■ 

h i 
e k i , where is the unit surface vector pointing towards 
cell k, and Am is the surface area connecting cells k and 
I. We then derive the following stochastic equations for 
mass and momentum exchange: 
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where dP k is the momentum exchange due to the fluctu- 
ating pressure tensor life, Vki = Vfe — v; and gkl is approx- 
imated on the surface kl by gki = \{pk + Pi)\{v k + v/). 
To close the discrete conservation equations we have to 
devise a discretization of the dissipative and fluctuat- 
ing parts which ensures the validity of the fluctuation- 
dissipation theorem. By choosing the discretization of 
the gradients d a 4> k — > A fc; e^0 fc /(2V fe ), the discrete 
momentum fluxes life and dPk take the form given in 
j 1 ll ] (see also supplementary information). The resulting 
set of stochastic differential equations Eqs. 1)112(1 may be 
integrated using various stochastic integration schemes 
|l5|: in this work we have used a simple Euler scheme. 

Molecular dynamics. The molecular descrip- 
tion is based on classical molecular dynamics and 
the CHARMM27 forcefield (incorporating the TIP3P 
paramctrization) which specifies bond, angle, dihe- 
dral and improper bonded interactions and non-bonded 
Lennard-Jones 6-12 and Coulomb interactions. The code 
is derived from a stripped down version of NAMD . 
We use a dissipative particle dynamics (DPD) thermo- 
stat [l7| ensuring local momentum conservation in such 
a way that hydrodynamic modes are not destroyed. 

Coupling protocol.- In our computational implementa- 
tion, the MD and FH components are independent cou- 
pled models || which exchange information after every 
fixed time interval At c . We set At c = ufh At = nMD^t, 
where At and St are the FH and MD time steps and, 
ufh and umd are integers which depend on the system 
being modeled; e.g. for water as solvent At c = 100 fs, 
n-FH = 10 and umd = 100. Conservation is based on 
the flux balance: both domains receive equal but oppo- 
site mass and momentum fluxes across the hybrid inter- 
face. This interface (H) uniquely defines the total sys- 
tem (MD+FH, see Fig. lb) and, importantly, the total 
quantities to be conserved. This contrasts with previ- 
ous schemes Q where particle and continuum domains 
intertwine within a larger overlapping region, preventing 
a clear definition of the system. 

The rate of momentum transferred across the hybrid 
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FIG. 1: The set-up used for our hybrid molecular simula- 
tions (a) and a close-up of the hybrid interface (b). The 
fluctuating hydrodynamics description (FH), resolved by the 
finite volume method, is coupled to a molecular model (MD) 
representing a dimyristoylphosphatidylcholine (DMPC) lipid 
monolayer solvated with water and restrained at the lipid 
head groups. We indicate by "P" and "C" respectively the 
particle and continuum cells adjacent to the hybrid interface 
"H". The buffer region of the MD system "B" (overlapping 
the C cell) is indicated by translucent water molecules and 
the water molecule density in the buffer region is shown in 
(c). 



interface is given by Fh = AJ H ■ e±_ , where ej_ is the unit 
vector perpendicular to the surface and the momentum 
flux tensor at "H" is approximated as J H = (Jp + J^)/2. 
Note that Jq involves the evaluation of the discretized 
velocity gradient at C, and thus requires the mass and 
momentum of the MD system at the neighbouring P 
cell averaged over the coupling time Atc- {Mp)& tc and 
(Pp)a* c , respectively (see Fig. lb). On the other hand, 
the momentum flux tensor at the P cell can be com- 
puted for the microstate using the kinetic theory formula 
J% = ([pVjVi + W l ]} A t c , with i 6 P and W% = J2j r u ' f y 
|18| being the contribution of atom i to the virial. Al- 
ternatively, Jp can be computed by introducing the 
coarse-grained variables at the neighboring MD and FH 
cells into the discretized Newtonian constitutive relation. 
Both approaches provide equivalent results in terms of 
mean and variance of the pressure tensor. 

The force Fh at the hybrid interface is imposed on 
the FH domain using standard von Neumann boundary 
conditions. In order to impose the force —Fh on the 
molecular system, we extend the MD domain to an ex- 
tra buffer cell ("B" in Fig. lb). Particles are free to 
cross the hybrid interface according to their local dy- 
namics, but any atom that enters in B will experience 
an external force —Fh/Nb which transfers the external 
pressure and stress. The number of solvent molecules at 
the buffer iVs(t) is controlled by a simple relaxation al- 
gorithm: AN B = ((N B )-N B )At c /TB, with t b ~ 500 fs. 
The average (Nb) is set so as to ensure that B always con- 
tains enough molecules to support the momentum trans- 
fer; here we use (Nb) = 0.75Mc/m, where Mc is the 
mass of the continuum cell C and m the molecular mass. 
Figure lc shows the equilibrium number density profile 
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of water at the buffer. Importantly, the density profile 
is flat around the hybrid interface. Due to the external 
pressure, it quickly vanishes near the open boundary. In 
fact, molecules eventually reaching this rarefied region 
in B are removed. If the relaxation equation requires 
ANb > 0, new molecules are placed in B with veloc- 
ities drawn from a Maxwellian distribution with mean 
equal to the velocity at the C cell. The insertion loca- 
tion is determined by the USHER algorithm [lfj, which 
efficiently finds new molecule configurations releasing an 
energy equal to the mean energy per molecule. Momen- 
tum exchange due to molecule insertion/removal is taken 
into account in the overall momentum balance . 

In fluid dynamics the mass flux is not an independent 
quantity but is controlled by the momentum flux [see 
Eqs. JU and (0)]. Consequently, we do not explicitly 
impose the mass flux on the MD system. Instead it arises 
naturally from the effect of the external pressure on the 
molecule dynamics near the interface. The mass flux J P H ■ 
ej_ is thus measured (via simple molecule count) from 
the amount of MD mass crossing the interface H over 
the coupling time At c . The opposite of this flux is then 
transfered to the adjacent C cell via a simple relaxation 
algorithm using a relaxation time (r r > 0(100) fs) 
large enough to preserve the correct mass distribution at 
the C cell, but still much faster than any hydrodynamic 
time. This guarantees mass conservation. 

Results. We first test the conservation of the total 
mass M and momentum P. Results are shown in Fig. 
2a, where we consider the equilibrium state of a hybrid 
MD simulation of water in a 3D periodic box 50 x 50 x 735 
A 3 (each cell is 50 x 50 x 15 A 3 ). The embedded TIP3P 
water domain (including the buffers) is 75A wide in the 
coupling (z) direction and was pre-equilibrated at 1 atm 
and 300K. Figure 2a shows the mean error in mass and 
momentum conservation. As stated above, mass conser- 
vation is ensured over a short time At c ~ O(100) fs, as 
clearly reflected in Fig. 2a. However, as the external 
force is imposed within the buffers B, the momentum 
conservation is ensured only on the "extended" system 
(MD+FH+B). The variation of momentum of the to- 
tal system (MD+FH) is then a small bounded quantity 
whose time average becomes smaller than the thermal 
noise after about lps (see Fig. 2a), i.e, faster than any 
hydrodynamic time scale. 

The FH description uses an accurate interpolated 
equation of state p(p) = (3.84 - 15.7p + 15. 3p 2 ) 10 4 bars, 
which fits for p = [0.54,0.70] g/mol/A 3 the outcome of 
NPT simulations of TIP3P water at T = 300K and pro- 
vides quasi-perfect match of the mean pressure, density 
(see Fig. 2b) and sound velocity. The shear and bulk 
viscosities of the FH model are assigned to match those 
of the MD fluid (for water at T = 300K we used those 
reported in Ref. Also, in cases where the viscosity 

varies locally, the FH model allows one to assign a dif- 
ferent viscosity for each cell. Momentum fluctuations at 




FIG. 2: (a) The normalized mean error in mass 
Em{t) and momentum Ep a (r) evaluated as E\(t) = 
(Ut° +T 5A(t)dt/T] 2 ) t0 /(M k } 2 , with 5A = A — (A). The 
dashed and solid horizontal lines are, respectively, the nor- 
malized standard deviation of mass and momentum within 
one cell (a[M k ]/ {M k ))- (b) Density field in a hybrid MD equi- 
librium simulation of water. Solid circles corresponds to MD 
cells. Error bars are the standard deviation of each cell den- 
sity. The grand-canonical (GC) result is (p) = 0.632g/mol/A 3 
and a[p] = 0.0045g/mol/A 3 . 



each cell are consistently controlled by the DPD thermo- 
stat in the MD region, and via the fluctuation-dissipation 
balance in the FH domain. Density fluctuations present a 
much more stringent test of thermodynamic consistency. 
Each fluid cell is an open subsystem so, at equilibrium, its 
mass fluctuation should be governed by the grand canon- 
ical (GC) prescription: a[p] = [pk B T/ (144)] 1 / 2 
(where a means standard deviation and c\ = [dP/dp)x 
is the squared sound velocity at constant temperature). 
Mass fluctuations within the MD and FH cells are both 
in agreement with the GC result (Fig. 2b ) in dicating 
that neither the usher molecule insertions [l0| nor the 
mass relaxation algorithm substantially alter the local 
equilibrium around the interface H. 

We now focus on transmission of sound waves which 
thus far have remained an open problem in the hybrid 
setting. In a slot of water between rigid walls we per- 
turb the equilibrium state with a Gaussian density per- 
turbation (amplitude 5% and standard deviation 45 A). 
As shown in Fig. 3a the resulting travelling waves cross 
the MD domain several times at the center of the slot. 
Sound waves require fast mass and momentum transfer 
as any significant imbalance would generate unphysical 
reflection at the hybrid interface. No trace of reflection is 
observed and comparison with full FH simulations shows 
statistically indistinguishable results. 

Finally, we validate the hybrid scheme against full MD 
simulations of complex fluid flow (set-up Fig. la). A 
sound wave generated by a similar Gaussian perturbation 
is now reflected against a lipid monolayer (DMPC) (Fig. 
3b) . Each lipid is tethered by the heavy atoms of the po- 
lar head group with an equilibrated grafting cross-section 
of 53 A /lipid, close to the experimental cross-section 
of membranes. In the hybrid simulation, the MD water 
layer close to the lipid membrane extends just 45A above 
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fluid rheology or slip flow past surfaces [3j. 
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FIG. 3: (a) Spatio-temporal diagram along z for the density 
field of a three-dimensional simulation of two sound waves 
traveling within a closed box filled with water. The region 
of width 45A around the centre of the box is described with 
molecular dynamics (MD), while the rest of the domain is 
solved via fluctuating hydrodynamics (FH). (b) The longi- 
tudinal velocity arising from the interaction between a sound 
wave and a grafted lipid layer (set up of Fig. la) . We compare 
hybrid MD (solid line), full MD (circles) and full FH simu- 
lation using purely reflecting walls (dashed line). Results are 
averaged over 15 nm from the monolayer; error bars indicate 
the standard deviation over 10 runs. 



it (see Fig. lb). Instead, in the MD simulation we consid- 
ered a large 180 x 50 x 50 A box of explicit water contain- 
ing around 50K atoms. The wave velocity near the layer 
is compared in Fig. 3b for the hybrid MD and MD simu- 
lations. The excellent agreement demonstrates that the 
coupling protocol accurately resolves features produced 
by the molecular structure. In Fig. 3b such effects are 
due to sound absorption by the lipid layer, highlighted 
by comparison with a FH simulation of the same wave 
impinging against a purely reflecting wall. The present 
sound waves simulations were done assuming an isother- 
mal environment. This is realistic if the rate of thermal 
relaxation Dxk 2 (with Dt ~ 1.5 x 10~ 7 m/s 2 the water 
thermal diffusivity and k = 2ir/X the wavenumber) is 
comparable with or faster than its sound frequency ck. 
The present simulations A ~ 50A are just in the limit 
of the isothermal sound regime |2ll |. while waves with 
A > O(10)A propagate adiabatically and require consid- 
eration of the energy flow |llj |. 

In summary, we have presented a stable and robust 
multiscalc method (hybrid MD) for the simulation of the 
liquid phase which embeds a small region, fully described 
by chemically accurate molecular dynamics, into a fluc- 
tuating hydrodynamics representation of the surrounding 
liquid. Mean values and fluctuations across the interface 
are consistent with hydrodynamics and thermodynamics. 
Sound waves propagating through the MD domain and 
flow behavior arising from the interaction with complex 
molecules are both treated correctly. We considered wa- 
ter waves reflected by DMPC monolayers, but the scope 
of this methodology is much broader, including inter alia 
the study of vibrational properties of hydrated proteins 
(via high frequency perturbatio ns) [j, y| , the ultrasound 
absorption of complex liquids j2f| or the simulation of 
quartz crystal oscillators for the study of complex 
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